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I. INTRODUCTION 

Due to its comparative simplicity, the short range bi- 
modal ±J Ising system is a widely studied model of a 
spin glass. The Hamiltonian is of the form proposed by 
Edwards and Andersoni 

H = - JijViVj (1) 

<ij> 

where the nearest-neighbor exchange interactions Jij are 
quenched random variables of fixed magnitude but ran- 
dom sign. These bonds are negative with a proba- 
bility p £ [0, 0.5] for the square lattice. The canon- 
ical model with p — 0.5 on the square lattice has 
been studied the most and it is well accepted that 
spin glass behavior occurs at zero temperature* 
and persists down to a critical probability p c of about 

q j ^8. 9. 10. 11. 12. 13. 14.15. 16. 17. 18. 19 

One of the compelling features of the 2-dimensional 
Ising model is the special property of allowing exact so- 
lutions - at least for finite systems in the absence of a 
magnetic field. A number of authors have taken this 
approach in various guiseo 13 i 17 i 20 i 21 ' 22 i 23 ' 24 i 25 i 26 . One ad- 
vantage of doing so is that it provides a direct access to 
the ground state properties without the need to extrapo- 
late from finite temperatures as is necessary, for example, 
in Monte Carlo based methods. 

Good agreement for the values of the ground state 
energy and entropy is obtained by the various workers. 
However, to the present, there has been very little devel- 
opment in extending the methodology to a direct calcu- 
lation of the spin correlations at zero temperature, and 
most of the results have come from extrapolations from 
finite temperatures. 

The present authors developed an approach 13,20,21 - 22 , 
based on the pfaffian matrix, that appears to capture the 
essence of the physics of the ± J system. The algorithm 
enables certain quantities such as the ground state free 
energy and entropy to be calculated exactly for very large 
lattices. The objective of the current paper is to extend 



that methodology to give direct access to the zero tem- 
perature correlation functions. 

Spin correlations for the bimodal Ising spin glass are 
expected to decay algebraically according to 

[< a a R > 2 ] av ~ iT" (2) 

at the critical temperature. In practice, finite size effects 
dictate that the spin correlations will decay exponentially 

[< <r a R > 2 ] av ~ JT" exp(-i?/£) (3) 

where the correlation length £ is expected to be propor- 
tional to the system size if the latter is large enough. 

To our knowledge the only direct computations of 
spin correlation functions in the ground state are those 
of Ozck: 27 . This work used a numerical transfer ma- 
trix method with long thin samples of circumference L 
wrapped around a cylinder of length 9L with open ends. 
The maximum possible circumference for this study was 
only L = 12, a consequence of the transfer matrix com- 
putational requirements scaling exponentially with L. 

All other attempts to study spin correlations in the 
ground state have involved extrapolation from finite tem- 
perature. Monte Carlo techniques have been employed to 
obtain results for low temperatures, for example 0.86 J 28 
and 0.63JS 9 .. However, it has never been clear just how 
reliable these extrapolations are. 

The transfer matrix method has also been used 
to obtain spin correlations functions at finite 
temperaturei^ii^ ., although the sample shape re- 
strictions are severe with cylindrical circumference 
L < 20. A better approach is probably the network 
modeli 7 - where the computational requirements scale 
as L 3 , not exponentially. Nevertheless, although larger 
values of L are feasible, the zero temperature limit is 
inaccessible. 

The algorithm we employ depends only on the number 
and distribution of frustrated plaquettes. This means 
that there is no need for the cylindrical circumference to 



2 



be especially small and the two lattice dimensions can 
be treated on an essentially equal footing. Although this 
algorithm is related in some formalities to the network 
model, it is especially designed to operate in the ground 
state. 

The method is also fully gauge invariant. This means 
that it is well suited for the determination of gauge in- 
variant quantities. In this context, a transformation of 
disorder is gauge invariant if the number and distribu- 
tion of frustrated plaquettes is unchanged. Examples of 
gauge invariant quantities are energy, entropy and the 
squares of correlation functions. In contrast, matching 
algorithms 24 - 25-26 , although more efficient, cannot deter- 
mine more than the energy. 

The formalism is developed in Section 2, and is then 
used in Section 3 for the evaluation of rj for the canonical 
± J model. 

II. FORMALISM 



where the matrix C is the same as D except for the re- 
ciprocal defects. 

The calculation of the partition function for the dis- 
ordered Ising model has been described before^. The 
key points are as follows. At zero temperature, D is a 
singular matrix with zero eigenvalues equal in number to 
the number of frustrated plaquettes. These eigenvalues 
(which occur in pairs) approach zero as some power of 
exp(-2J/fcT) 

e = ±~X cxp(-2Jr/kT) (7) 

where r is an integer. The ground state energy and en- 
tropy can be expressed exactly as 

F = -2J + 2J^r d (8) 

d 



The planar Ising model has long been analytically ac- 
cessible since it can be mapped onto a system of non- 
interacting fermions. This mapping can take various 
forms, including the transfer matrix and combinatorial 
methods. The combinatorial, or Pfaffian, method^ is 
particularly well suited to the study of disordered planar 
Ising systems. Essentially, each lattice site is decorated 
with four fermions which have interactions across bonds 
as well as intrasite interactions. 

The partition function for a disordered planar Ising 
model can be expressed in the form 



Z = 2 



iV 



[ JJ cosh(Jy/AT)] (detZ)) 1 / 2 



(4) 



<ij> 



where the product is over all nearest neighbor bonds Jjj 
on the N site lattice and D is a skew-symmetric matrix. 
The square root of the determinant of D is the major 
feature and is precisely the Pfaffian^i. It is also propor- 
tional to the trace over all closed lattice polygons, that 
is 



Tr Y[ (1 + tijaiaj) = 2- N (detD) 



1/2 



(5) 



<1J> 



where ty = tanh(Jy/feT). 

The correlation functions can be expressed within the 
same formalism as a reciprocal defect proble m 2 - ' -3.1 i . We 



choose a path between a pair of spins and replace Uj with 
tZ- for bonds on that path. It may be noted here that 
the path can also be disjoint, in which case the formalism 
would give a correlation function for four or more spins. 
In terms of determinants, the correlation function can be 
easily expressed as 



< <j vr > 



detC TT t 2 
detD 11 ij 

path 



(6) 



S = k^2lnX d 



(9) 



The sums are over all zero eigenstates. An algo- 
rithm based on degenerate state perturbation theory was 
developed^ to evaluate exactly the r and X, and hence 
the ground state energy and entropy. In this paper we 
show how to extend this algorithm to the correlation 
functions as well. 

It was found^i to be convenient to transform the lattice 
fermions from the sites to the bonds. In this way, we 
can associate two fermions with each bond and, on the 
square lattice, four fermions with each plaquette. Figure 
^ provides an illustration. A simple generalization to 
another planar lattice has proved straightforward^. 
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FIG. 1: A plaquette showing how lattice fermions are asso- 
ciated with bonds and plaquettes. 

For the ± J model in particular we write 



D = Dn 



5D l 



(10) 



where S = 1 — t with t = tanh(J/fcT). This equation 
is exact and D\ characterizes a perturbation away from 
the singular matrix Dq and has nonzero matrix elements 
only across bonds. The singularities of -Do arise due to 
the frustration. For each frustrated plaquette, Dq has 
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a zero eigenvalue with an eigenvector localized on the 
four associated fermions. To determine the ground state 
properties, it is necessary to determine the defect eigen- 
values of D, that is those eigenvalues that are zero at zero 
temperature. First D\ is diagonalized in the basis of the 
defect eigenvectors localized on the frustrated plaquettes. 
The second order calculation then involves diagonalizing 

D 2 = D igc D x (11) 

in the basis of the eigenvectors corresponding to zero 
eigenvalues at first order. This process is continued order 
by order where, for n > 2, 

D n = D n _i(l + G n _ 2 D n _ 2 ) •••(! + G x D t )g c D x (12) 

until no zero eigenvalues remain. In these equations, g c 
is the continuum propagator and G n is the propagator 
for eigenstates determined at order r»21. 
The matrix C is given by 

C = D + 5(D 1 +2V) + (6 2 + 6 3 + ---)V (13) 

where V = — D\ for matrix elements across a path bond. 
All other elements of V are zero. This equation is easily 
derived by expanding t^ 1 —t in powers of S. Since t 2 = 1 
in the ground state, we can now state our main task as 
the computation of the limit 

r 2i i- detC . . 

< <to°r > t=o = Jim 14 ) 
<5^o det D 

To achieve this goal efficiently, we have discovered a 
remarkable property of the continuum propagator. It 
can be written as the sum of two terms 

9c = 9d + 9c2 (15) 

where g c \ is 4 x 4 block diagonal in the four fermions as- 
sociated within a plaquette and g c2 is 2 x 2 block diagonal 
in the two fermions associated with a bond. The physical 
relevance of this decomposition can be realized in that g c2 
should not play any role in the ground state, being only 
an issue for excited states. Both g c2 and D\ are bond 
diagonal, meaning that g c2 D\ costs energy while taking 
us nowhere. We can now state an important theorem 
concerning the determination of det C . 

Theorem. In the ground state limit, the solution for 
det C using 

C = D + S{D 1 +2V) + (5 2 +S 3 + ■■■ )V 

g c = g c i + g C 2 (16) 

is exactly the same as the solution using 

C = D + 6(Di + 2V) 

9c = 9d (17) 



Essentially, the effects of g c2 exactly cancel the awk- 
ward nonlinear terms in 5. We can simply rerun the 
perturbation theory for D\ + 2V instead of D\, noting 
that this simply requires a sign change across path bonds. 
Further, there is an obvious corollary that says that g c2 
can be simply ignored while computing det_D. The con- 
tinuum propagator is really just g c \ and is localized inside 
plaquettes. The physical sense of this is clear. A proof 
of the theorem is given in the appendix. 

In2l it was proven that det D is proportional to the 
ground state degeneracy (Equation © is an equivalent 
statement). We can generalize this by considering detC 
as proportional to an effective degeneracy for the recip- 
rocal defect system. In fact we can usually write 

<ct ct r > 2 =cM2(Sc{R)-S d )) (18) 

where Sd is the system entropy and Sc is the analog for 
the reciprocal defect system. This result applies unless 
the effective energy of the reciprocal defect system is less 
than the actual energy in which case the correlation func- 
tion is zero. We are not aware of this result appearing 
elsewhere in the literature. As a function of system size 
L, the entropy is expected to vary as^ 

S D (L) = S m +S D2 /L (19) 

with Sd2 > for p > p c and Sdi < otherwise. This 
same behavior is also expected with cylindrical winding, 
that is periodic boundary conditions applied in one di- 
mension. From this, we can of course reasonably expect 
that the first cumulant approximation 

< o~qo~ > 2 ~ A{R) exp(-B(R)/L) (20) 

is valid if L is sufficiently large. This is of course in accor- 
dance with finite size scaling theory, where it is expected 
that A(R) - RT^ and B(R) ~ R if R is sufficiently large. 
The nature of higher order corrections with cylindrical 
winding is unclear—. 

III. RESULTS 

We have first estimated correlation functions with 
paths of length R, in units of the lattice spacing, parallel 
to the axes of long thin cylinders of dimensions 9L x L 
with L — 12, 16, 32 and 64. OzekiSi used similar ge- 
ometry, but the largest size treated was L = 12. The 
results are displayed in figure El The error bars indicate 
uncertainties equal to two standard deviations, that is a 
95.4 per cent confidence interval. The averages are over 
10 5 random samples, except for L = 64 where only 10 4 
samples were obtained. The uncertainties do not depend 
strongly on either L or R, only on the number of sam- 
ples. With 10 5 samples, they all sit in the range 0.0021 



to 0.0026. The corresponding interval with 10 4 samples 
is 0.0073 to 0.0086, about a factor of /L0 larger. The 
curves are fits of the data to the form 



< a a R > 2 = AFT 11 exp(-R/g) 



(21; 
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These \ 2 nonlinear fits were done using the Levenberg- V 
fyfeirquardt incthor 133 with the statistical uncertainties i,n 
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FIG. 2: The square of the spin corrBlation function on 9L x V" 
lattices with L = 12 (pluses), 16 (crosses), 32 (stars) and 64 
(squares). 

The coefficient A and the exponent ij are both approx- 
imately independent of L. The values of A are 0.73(4), 
0.72(2), 0.72(1) and 0.71(3) respectively for L = 12, 16, 
32 and 64. The corresponding values of r] are 0.12(6), 
0.14(3), 0.15(1) and 0.14(2). In comparison, the correla- 
tion length £ does not vary like L as would be expected 
from finite size scaling theory. In fact we find that it 
varies more like . The respective values of are 
0.41(6), 0.41(4), 0.40(3) and 0.47(10). We also find that 
finite size corrections to the entropy for the 9LxL system 
scale like L~?, whereas scaling as equation (|19|1 applies 
with the L x L lattice. Clearly the correlation functions 
are intimately associated with the entropy as is high- 
lighted in equation i|18|) . Interestingly Lukic et al2i also 
observe finite size corrections to the ground state entropy 
scaling like L~2 ( albeit for a system with different geom- 
etry. ^ 

Figure shows the same data plotted against L 2 for 
R = 4, 8, 16 and 32. The x 2 fits are to the function 
a(R) exp(— P(R)L~z ) and were first done with linearized 
data before polishing with nonlinear fits to obtain uncer- 
tainties for the fitting parameters a(R) and (3(R). A 
consequent fit of a(R) to a power function R _r! gives 
i] = 0.14(1). We have also tried to fit the data to the 
form of equation (|20l) . The result is certainly inferior 
with considerably larger \ 2 values for R = 8 and 16 in 
particular, although B(R) does follow R rather roughly. 
In detail, for R = 4, 8, 16 and 32, the respective values 
of B(R)/R are 0.65(6), 0.69(4), 0.72(4) and 0.65(9). 
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FIG. 3: As figure |3 with the data plotted against for 

4 (pluses,), 8 (crosses), 10 (stars) and 32 (squaies). 
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FIG. 4: The square of the spin cofrelation function on L x 
L lattices with L — 32 (pluses), 48 (crosses), 64 (stars), 96 
(squares) and 128 (filled squares). 



We have also performed estimations of correlation 
functions on L x L lattices. The lattices are cylindri- 
cally wound as before but the reciprocal defect path is 
around the circumference of the cylinder, as far as pos- 
sible from the nested boundaries. Figure 01 shows the 
results for a range of values of L up to 128. The number 
of random samples taken was 10 5 . The curves are a fit to 
the form of equation (|21|) and serve only as a guide to the 
eye. The statistical uncertainties for our estimates, with 
10 5 random samples, on L x L lattices are in the interval 
0.0021 to 0.0026, the same as for the 9L x L lattices. 

Figure [S] shows data plotted against L~ l , The fits are 
to the form of equation l|20|l . Only data for which L > 
2R is used since it is observed that data for smaller L 
sits significantly above the fit. This is most probably 
a consequence of the finite size effect of the cylindrical 
circumference. A further fit of A(R) to R~ v reveals 77 = 
0.13(1) which is in reasonable agreement with the values 
found above. 
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FIG. 5: The square of the spin correlation function on L x L 
lattices with R = 4 (pluses), 8 (crosses), 16 (stars) and 32 
(squares). 



The eigenvalue equation is formally written 



C | *) = A | *) (A2) 
A and | W) can be expanded in powers of 5 



X = S n ei + 0(5 n+L ) 



*) =i *o) + <n *o 



(A3) 



(A4) 



where the eigenvalue is the jth defect eigenvalue of order 
with n > and | Wo) =| n,j) is the leading term for 
thd corresponding eigenvector, that is D \ n,j) = 0. 
The notation here is in keeping withSi. 

For 1 < m < n, equating powers of 5, we obtain 



IV. DISCUSSION 

Earlier work on the exact calculation of the energy 
and entropy of the 2-dimensional ± J spin glass has been 
extended to cover spin correlations in the ground state. 
The computational requirements of the algorithm depend 
only on the number and distribution of frustrated plaque- 
ttes, and this provides the possibility of studying quite 
large samples. The big advantage of the approach is that 
there is no need to extrapolate from finite temperatures 
as in Monte Carlo methods. 

Our best estimate of the exponent is rj — 0.14. We 
regard the cylindrical winding as the most reliable con- 
figuration for calculation since it provides one long di- 
mension. This value for rj emerges from the calculation 
for the largest value of L, and for the composite analysis 
in Figure |3| The calculation on the L x L lattice also 
yields a similar figure. 

A wide range of values of rj have been reported in the 
literature over the years. A significant number of these 
favor a value around 0.2, somewhat higher than ours, but 
most of these use extrapolated Monte Carlo data rather 
than focusing directly on the ground state. Interestingly 
a recent calculation 7 reports ry = 0.138 ± 0.005 in agree- 
ment with the current calculation. 



APPENDIX A 

Although the central theorem of this article makes 
good physical sense, it can also be proven with math- 
ematical rigor. This appendix outlines a proof. For 
present purposes, all matrices are defined as pure imagi- 
nary Hermitian, that is with real eigenvalues. 

The main issue is finding the eigenvalues of 



C = D a + S(D 1 + 2V) + (S 2 + S 3 + ■ ■ ■ )V (Al) 



p=2 

(A5) 

where C\ = D\ + 2V and the sum vanishes for m < 2. 
^From this it follows that 



(r,i | d | * m _ x ) + ^2(r,i \ V | <F m _ p ) = 8 mn 5 rn 8 l0 e j n 

p=2 

Then, for m = n = r = 1 we now have 



|d | l,j) = <%ej 



(A7) 



and the first order states are determined by diagonaliza- 
tion of d in the defect basis set, that is the localized 
vectors associated with the frustrated plaquettes. 

Following a development similar to that given for the 
matrix D inSi we can show that 



*m) = ffcd I *m-l> +9cV^2 I *m-p> 

p=2 

S m ng c ei \n,j) + 22\ r, i)Z l r (A8) 



where Z l r is an undetermined coefficient. 

For 1 < m < n — 1, we now use equations ()A6|) and 
(IA8I) to proceed to the second order problem. As stated 
in the main text, we write g c = g c \ + g C 2- Further, it 
is easy to demonstrate that g C 2D\ — |, Vg C 2V = — hV, 
Cig c2 V = -\V and C 1 g c2 C 1 + V = \C X . Then, defining 
C*2 = dffcidj it can be shown that, for r > 1, 



- m—p/ 



r,i\C 2 | W m _i) +^>,z | C x g cX V \ W, 

p=2 

= &m,n-l$Tnkj4i (A9) 
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and we can determine the second order states using 



(2,i\C 2 \2,j)=5„ 



*J C 2 



(A10) 



Also, with r = 1, we can now find the coefficients 



1 - 

(All) 



p=2 



and. with the definition 



(A12) 



we find that 



I *m) = ((l + ^lCi) ffc i+.9 c2 ) 



Cl I *m-l> + ^E I 
V P=2 / r>l,t 



E i r > j >^ 



(A13) 



Note that we are using F r for the propagators of the 
system with reciprocal defects to avoid confusion with 
the G r used in2i for the diagonalisation of D. 

Next we consider 1 < m < n — 2 and work with equa- 
tions 1|A9|) and ilA13(l as well as the results C 2 g c2 V = 
— \C\g c \V and C 2 g c2 C\ — ^C 2 — Cig c \V. We also de- 
fine C 3 = C 2 (l + FiCi)g c iCi. Then, for r > 2, we can 
arrive at 



(r,t | C fe | * m _i)+ 
£)<r, i | C fc _i(l + F fc _ 2 C fe _ 2 ) • • • (1 + FiCi)ff c iV | * m -p) 

p=2 

= 5m,n-k+l5rn&ije! n (A17) 

and, second 



| * TO ) = ((1 + F fc _iC fc _i) ■ • • (1 + FiCiJffd + <? c2 ) 

- E i 



m—p 

V p=2 



(A18) 



We can note that, comparing with equations l|A14|) and 
HA15jl . these assumptions are both true for /c = 3. Also, 
with the definition (|A16|) . we can easily show that 



C k g c2 V = 
and 



--C fc _i(l + F fe - 2 C fe _ 2 ) • • • (1 + FiCO^y 

(A19) 



Cfeff^Ci = ^C fc -C fc _i(l+F fc _ 2 C fc _2) • • • (l+FiCiJffdV 

(A20) 

Now, for 1 < m < n — k and r > fc, we can use the 
assumptions (|A17|) and (|A18fl with 1A19|) and (|A20|) to 

prove two results. First, for r > k, 



(r,i | C 3 | E(M I C 2 (l + FiCi) 5 ciV I *m- P ) 

p=2 

=^m,Ti-2^rn^e^ (A14) 

Clearly the third order problem is solved by diagonalizing 
C3. Furthermore, for r = 2, it is possible to determine 
the coefficients Z\ and we find that 



I * m > = ((1 + F 2 C 2 )(1 + *i<?l)ffcl + 5c2 ) 

x j Ci I *m-i) + vjr I * m _ p ) j + E I r^Z* 

V p=2 / r> 2 ,i 

(A15) 

To prove the main result, we can use induction. First 
we define, for k > 3, 



C k = C fc _i(l + F fe _ 2 C fc _ 2 ) • • • (1 + fiCO^d (A16) 

and, for 1 < m < n — k + 1 and r > fc — 1, introduce two 
assumptions. First, 



(r,i I C k+1 |* m _i/ 



X>,<| C fc (l + ^_ 1 C , fe _ 1 )---(l + J F 1 C 1 ) ffcl 1/ I * m _ p ) 

p=2 

= Sm^n-kSrnSije^ (A21) 

and, second, after determining the coefficients Z l k , 
I * m > = ((1 + F k C k ) •••(! + FiCiJfld + <? c2 ) 



d 1 M+^E 1 

V P=2 

E i r '^ 



(A22) 



r> k ,i 



This completes the proof of the theorem. Clearly the 
assumptions (|A17|I and l|A18l) imply the results (|A21|) 
and l| A22|l . The eigenvalues of C are determined by di- 
agonalisations of the C r which are defined independently 
of g c2 - 
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